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ABSTRACT 


Three-dimensional  models  of  the  crust  and  mantle  structure  beneath  Area  1  International  Monitoring 
System  stations  in  Eurasia  are  constructed  for  use  with  either  asymptotic  ray  or  numerical  methods  of 
waveform  modeling.  The  models  combine  crustal  and  upper  models  of  varying  resolution  specified  on 
latitude  and  longitude  grids  having  variable  spacing.  Model  parameters  are  interpolated  using  Delaunay 
triangulation  in  3-D  (tetrahedra).  Unless  accurate  narrow-angle  crustal  reflections  and  reverberations  are 
required,  ray  bookkeeping  is  simplified  by  making  all  first-order  crustal  discontinuities  narrow  transition 
zones.  Station-specific  path  corrections  (SSPCOs)  for  travel  times  are  obtained  in  these  models  by  dynamic 
ray  tracing  (DRT).  DRT  provides  information  on  wavefront  that  can  be  used  to  accurately  interpolate  travel 
times  by  a  paraxial  approximation  in  the  vicinity  of  end  points  of  rays.  By  this  method  travel  times  are 
computed  in  3-D  models  at  dense  grids  specified  around  each  IMS  station.  Each  grid  point  will  contain 
travel  time  and  quantities  needed  to  interpolate  travel  times  spatially  at  finer  intervals.  These  interpolating 
quantities  can  also  be  useful  to  event  relocation.  The  grid  of  values  is  given  as  a  binary,  direct  access  file. 
The  wavefront  curvature  information  provided  for  each  path  includes  information  needed  for  path 
integrated  attenuation  (t*),  geometric  spreading,  and  ray-synthetic  seismograms. 
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OBJECTIVE 


Travel-time  tables  will  be  generated  for  paths  of  regional  seismic  phases  within  and  to  area  1  IMS  stations  in 
Eurasia  using  three-dimensional  models  of  crustal  structure.  The  task  assigned  to  the  University  of  Connecticut 
is  to  develop  algorithms  for  constructing  three-dimensional  models  of  the  crust  and  upper  mantle  and  to 
compute  travel  times  and  ray  theoretical  waveforms  for  regional  seismic  waves  propagating  in  these  models. 
The  results  of  this  work  will  provide  accurate  regional  phase  travel  times  for  improvements  in  the  IDC’s 
location  estimates  based  on  data  from  IMS  stations  in  Eastern  Asia.  The  incorporation  of  three-dimensional 
structure  will  correct  for  the  effects  of  off-azimuth  paths  on  travel  times,  which  are  pervasive  throughout  this 
region. 


RESEARCH  ACCOMPLISHED 


Method  of  computing  travel  times  and  SSPC’s 

Under  work  completed  by  the  co-Principal  Investigator  under  an  AFOSR  contract,  a  three-dimensional  dynamic 
ray  tracing  program  was  written  to  enable  prediction  of  Lg  from  SmS  ray  paths  in  three-dimensionally  varying 
crustal  models.  Rays  were  traced  by  integrating  kinematic  and  dynamic  ray  tracing  equations  (Cerveny  and 
Hron,  1980;  Cerveny,  1985).  The  dynamic  ray  tracing  system  gives  information  needed  to  compute  wavefront 
curvature,  which  can  be  used  to  compute  geometric  spreading  and  quantities  needed  for  summation  of  Gaussian 
beams,  a  technique  of  seismogram  synthesis  closely  related  to  the  Maslov- WKBJ  method.  These  techniques 
have  been  applied  by  the  co-PI  to  a  wide  variety  of  applications  requiring  two-point  ray  tracing  (e.g.,  Cormier 
and  Beroza,  1987)  and  ray  theoretical  synthesis  of  seismograms  in  three-dimensionally  varying  media  (e.g., 
Cormier  and  Anderson,  1999;  Cormier,  1986). 

Three-dimensional  dynamic  ray  tracing 

Station-specific  path  corrections  (SSPC’s)  will  be  obtained  from  travel  times  computed  by  dynamic  ray  tracing 
(DRT).  DRT  provides  information  on  wavefront  that  can  be  used  to  accurately  interpolate  travel  times  by  a 
paraxial  approximation  in  the  vicinity  of  end  pints  of  rays.  By  this  method  travel  times  will  be  computed  in  3-D 
models  at  dense  grids  specified  around  each  IMS  station.  Each  grid  point  will  contain  a  travel  time  and 
quantities  needed  to  interpolate  travel  times  spatially  at  finer  intervals.  These  interpolating  quantities  can  also 
be  useful  to  event  relocation.  The  grid  of  values  will  be  given  as  a  binary,  direct  access  file,  and  will  be  the  first 
deliverable  product  of  the  work.  The  wavefront  curvature  information  provided  for  each  path  will  also  provide 
information  needed  for  path  integrated  attenuation  (t*),  geometric  spreading,  and  ray-synthetic  seismograms. 

Travel  times  for  arbitrary  3-D  paths 

The  second  stage  in  this  effort  will  consist  of  a  simple  code  that  retrieves  the  ray  end-point  quantities  from  that 
file  and  calculates  the  travel  time  of  a  specified  phase  from  a  source  by  a  paraxial  approximation.  An  important 
advantage  of  the  paraxial  approximation  is  the  avoidance  of  expensive  two-point  ray  tracing,  which  must  find 
the  exact  ray-path  connecting  source  and  receiver  to  high  accuracy.  This  will  enable  rapid  retrieval  of  travel 
times  for  any  arbitrary  path  in  a  fully  three-dimensional  structure,  including  the  effect  of  three-dimensional  path 
deviation  from  a  great  circle.  In  the  theory  of  dynamic  ray  tracing,  the  paraxial  approximation  calculates  the 
travel  time  from  position  xG  to  position  x  by 


T(x,xJ  =  T(x)  +  p-Ax  +  H  Ax  MAx' H1 


where  p  is  a  vector  slowness  at  the  ray  end  point,  M  is  a  3  x  3  matrix  of  second  derivatives  of  the  travel  time 
field  in  ray-centered  coordinates  at  the  ray  end  point,  x  is  the  Cartesian  vector  difference  between  the  station 
location  within  a  grid  box  and  the  ray  end  point  in  the  grid  box,  and  H  is  a  3  x  3  matrix  whose  columns  consist 
of  the  vectors  of  the  ray-centered  coordinate  basis.  The  algorithm  of  the  proposed  code  simply  evaluates  the 
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above  equation  using  25  quantities  stored  in  grid  boxes  for  each  IMS  station  corresponding  to  paths  from  all 
other  boxes  not  containing  that  station.  It  is  simple  to  include  among  these  quantities  the  path  integrated 
attenuation  (t*)  for  regional  phases  using  measurements  from  models  determined  from  the  DSS  data  analyzed 
by  the  University  of  Wyoming  group. 

Ray  paths  need  not  precisely  connect  source  and  receiver.  The  paraxial  equation  above  accurately  corrects  the 
travel  time  to  the  time  that  would  be  computed  for  the  exact  ray,  provided  that  the  source  and  receiver  are  not 
too  distant  from  a  frequency-dependent  region  of  validity  often  termed  the  “paraxial  vicinity.”  In  tests  of  the 
algorithm  for  the  models  and  ranges  proposed  here,  use  of  the  paraxial  approximation  at  ray  end  points  to 
extrapolate  travel  times  within  25  km  x  25  km  boxes  were  found  to  produce  travel  times  accurate  to  within  0.01 
sec. 

Previous  work  (Cormier  and  Anderson,  1999)  modeled  the  Moho  and  the  boundary  between  a  sedimentary 
layer  and  hard-rock  layer  (basin  topography)  as  surfaces  interpolated  by  splines  under  tension.  The  Moho  and 
basement  surfaces  were  assumed  to  spatially  continuous  transition  zones.  With  surfaces  of  discontinuities  so 
modeled,  all  rays  become  either  (1)  turning  rays  and  multiply  surface  reflected  rays  within  the  model  or  (2)  rays 
that  dive  deeply,  leaving  the  bottom  of  the  model.  This  procedure  eliminates  the  need  for  computing  many 
reflection  coefficients  and  complex  descriptions  of  rays  interacting  with  multiple  first-order  crustal 
discontinuities. 

Starting  three-dimensional  models 

Starting  models  for  three-dimensional  structure  were  assembled  from  published  models  for  the  crust  and  upper 
mantle  of  Eastern  Asia  (Fielding  et  al,  1992;  Barazangi  et  al.,  1996;  Mooney  et  al,  1998;  Ritzwoller  and 
Levshin,  1998;  Sambridge  and  Gudmundsson,  1998).  Each  of  these  models  were  constructed  to  satisfy  different 
types  of  data  and  reported  with  differing  degrees  of  resolution.  Hence,  there  is  no  guarantee  at  the  outset  that  a 
hybrid  model  incorporating  features  of  all  these  models  will  be  successful  in  reducing  the  travel  time  residuals 
from  ground  truth  events  calculated  from  laterally  homogeneous  reference  models. 

The  Cornell  model  of  Barazangi  et  al.  (1998)  reports  the  depth  to  the  Moho  and  depth  to  a  hard-rock  basement 
in  the  crust  at  0.1  x  0.1  deg  resolution.  Crust  5.1  by  Mooney  et  al.  (1998)  incorporates  many  published  local 
and  regional  refraction  studies  in  a  global  7  layered  crust  model  reported  at  either  a  1  x  1  deg  resolution  for 
crustal  and  basement  thickness  or  a  5  x  5  deg  resolution  for  the  detailed  7  layered  model.  The  University  of 
Colorado  Model  (or  CU  model)  of  Ritzwoller  and  Levshin  (1998)  is  obtained  from  surface  wave  group 
velocities  and  some  Pn  velocities.  Except  for  the  P  velocity  at  the  underside  of  the  Moho  constrained  by  Pn,  the 
P/S  velocity  ratio  is  not  well  constrained.  Although  Sambridge  and  Gudmundsson’s  (1998)  Regionalized  Upper 
Mantle  (RUM)  model  is  based  on  a  tectonic  age  regionalization  (Jordan,  1981)  at  a  spatially  fine  scale,  it  has 
been  inverted  from  primarily  teleseismic  data.  RUM,  as  it  name  implies,  is  a  model  of  the  upper  mantle  rather 
than  of  the  crust. 

Since  previous  work  by  Cormier  and  Anderson  (1999)  demonstrated  the  importance  of  fine  scale  moho  (10  km 
and  less)  basement  topography  on  the  progagation  of  Pn  and  Lg,  it  was  decided  to  incorporate  into  a  hybrid 
model  the  high  resolution  model  (0. 1  deg  x  0. 1  deg)  of  basin  thickness  and  Moho  topography  of  the  Cornell 
model.  Since  Crust  5.1  provides  high  detail  in  the  vertical  direction  (up  to  7  layers),  it  was  used  for  structure 
above  the  Cornell  high  resolution  Moho.  (The  high  resolution  basin  structure  of  the  Cornell  model  has  yet  to  be 
included  in  test  models.)  The  RUM  model  has  been  used  for  structure  below  the  Cornell  Moho. 

Parameterization  and  interpolation  of  three-dimensional  models 

The  varying  types  of  data  for  crust  and  upper  mantle  structure,  collected  at  widely  different  spatial  scales,  and 
the  highly  uneven  distribution  of  ground  truth  events  present  a  challenge  to  the  parameterization  a  three- 
dimensional  model  for  travel  time  computation.  The  chosen  model  parameterization  should  be  flexible  enough 
to  be  specified  at  high  resolution  where  data  is  available  and  at  lower  resolution  where  it  is  not.  Resolution 
should  be  high  to  describe  features  important  to  regional  wave  propagation,  such  as  Moho  and  basin 
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topography,  but  can  be  lower  near  interfaces  having  smaller  velocity  contrasts  and  lower  with  increasing  depth 
in  the  mantle,  where  heterogeneity  power  decreases.  To  allow  a  smooth  transition  from  the  regional  times  and 
waveforms  to  teleseismic  times  and  waveforms,  the  parameterization  should  incorporate  the  near  sphericity  of 
earth  and  its  major  structural  discontinuities.  The  parameterization  should  also  be  economical  in  the 
specification  of  the  number  of  knots. 

The  proposed  algorithm  of  kinematic  and  dynamic  ray  tracing  requires  known  first  and  second  spatial 
derivatives  of  a  3-D  model,  which  can  be  met  by  the  use  of  spline  interpolation  using  polynomials  of  order  3  or 
higher.  Since  the  investigator  had  previously  applied  splines  to  interpolate  model  parameters  in  3-D,  tests  were 
performed  to  determine  whether  this  type  of  interpolation  could  also  be  applied  to  the  3-D  structure  beneath 
IMS  stations.  The  routine  previously  used  with  dynamic  ray  tracing  employed  splines  under  tension  at  equally 
spaced  knots  in  Cartesian  co-ordinates  (Cormier,  1986).  The  FITGRID1  code  previously  employed  under 
tension  code  was  not  easily  adaptable  to  our  IMS  model.  Its  restriction  to  equally  spaced  knots  requires 
cumbersome  and  constant  resampling  to  incorporate  both  earth  sphericity  and  spatially  variable  resolution.  In 
future  experiments,  the  software  described  by  Wessel  and  Bercovici  (1998),  which  may  allow  for  irregularly 
spaced  3-D  data,  will  be  tested  for  application  to  dynamic  ray  tracing. 

The  most  recent  release  of  the  NCAR  graphics  package2  has  a  set  of  routines  named  CSAGRID3  that  can 
interpolate  by  splines  on  an  irregularly  spaced  Cartesian  grid.  Unfortunately,  tests  showed  that  without  an 
available  tensioning  factor  in  the  current  version  of  CSAGRID,  the  splines  very  poorly  tracked  realistic  seismic 
velocity  variations  within  layers.  This  deficiency  could  be  remedied  only  be  expensive  densification  of  knot 
points,  ultimately  making  its  use  as  cumbersome  as  spline  routines  having  regularly  spaced  grids. 

The  next  and  preferred  parameterization  and  interpolation  scheme  tested  is  the  one  used  and  advocated  by 
Sambridge  et  al.  (1995)  in  the  construction  of  the  RUM  model.  This  scheme  parameterizes  3-D  earth  models  by 
knots  connected  by  tetrahedra.  A  linear  gradient  in  velocity  is  assumed  for  the  interpolated  quantity  within  each 
tetrahedral  element,  making  it  possible  to  analytically  integrate  within  each  tetrahedron  both  the  kinematic  and 
dynamic  (geometric  spreading  and  wavefront  curvature)  ray  tracing  equations.  A  public-domain  software 
package  named  qhull  (Barber  et  al,  1996),  maintained  by  the  Computational  Geometry  Center  at  the  University 
of  Minnesota4  ,  is  available  for  performing  the  Delaunay  tetrahedral  tesselation  of  a  model  volume  given  a  set 
of  knots  and  model  values.  The  file  containing  model  values  at  knots  may  be  given  with  knots  occuring  at  any 
arbitrary  order  or  spacing.  Qhull  returns  a  list  of  index  numbers  of  the  4  knots  describing  each  tetrahedron. 
Routines  for  navigating  though  the  tetrahedra  are  available  from  Sambridge  5 

Figure  1  shows  an  example  of  a  hybrid  3-D  model  constructed  from  qhull  tetrahedrall  tesselation  using  the 
Cornell  Moho,  the  Crust  5.1  crust,  and  the  RUM  mantle.  The  example  is  for  a  200  x  200  x  200  km  block  whose 
surface  is  centered  at  Nilore,  Packistan.  Note  the  thickening  of  the  Moho  to  the  north  of  Nilore,  an  effect  of  the 
crustal  root  that  has  grown  from  Indian-Eurasian  plate  collision.  Some  spatial  variability  is  also  seen  in  the  low- 
velocity  zone  given  by  RUM.  The  sampling  of  the  low  velocities  in  the  upper  most  crust  2.5-4km/sec  needs  to 
be  densified  to  remove  some  artifacts  associated  with  too  broad  tetrhedron  facets  at  the  surface  of  the  earth. 

Ray  tracing 

Although  the  tetrahedral  tesselation  of  the  model  allows  analytic  integration  of  the  ray  tracing  equations,  intitial 
tests  of  travel  times  have  used  a  Runge-Kutta  numerical  integration  to  shoot  rays.  Figure  2.  Shows  an  example 
of  ray  end  points  shot  at  uniform  increments  vertical  and  azimuthal  take-off  angles  from  NIL.  Note  some 
effects  of  three-dimensional  structure  are  evident  in  the  scatter  of  rays  at  longer  distances.  A  strong  shadow 
zone  is  seen  around  10-15  degrees.  In  this  distance  range,  the  first  arrival  is  a  Pn  phase  just  grazing  the 
underside  of  the  Moho. 

CONCLUSIONS  AND  RECOMMENDATIONS 


Additional  hybrid  3-D  models  will  will  be  constructed  for  testing  SSPC’s  for  IMS  using  the  a  tetrahedral 
tesselation  described  in  this  report.  A  denser  sampling  will  be  incorporated  in  the  upper  2-3  km  of  the  crust  to 
better  describe  the  known  three-dimensional  variations  in  soft  sediment  cover.  To  speed-up  forward  modeling 
of  travel  times  and  dynamic  ray  tracing  quantities,  an  analytic  integration  (Menke  and  West,  personal 
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communcations)  for  kinematic  and,  eventually  dynamic  ray  tracing  within  tetrahedra,  will  be  substituted  for  the 
numerical  integration  currently  used.  A  test  of  the  spline  under  tension  parameterization  of  Wessel  and 
Bercovici  (1998)  will  also  made  to  see  if  it  offers  any  greater  flexibility  in  model  parameterization  compared  to 
the  tetrahedral  parameterization. 

It  will  also  be  important  to  test  the  accuracy  of  travel  time  calculations,  which  may  be  affected  by  the  choice  of 
parameterization  and  the  accuracy  of  paraxial  approximations  against  other  methods.  A  test  will  be  performed 
with  at  least  one  model  comparing  the  results  of  dynamic  ray  tracing  and  the  results  of  the  numerical  eikonal 
solution  method  of  Vidale  (1988)  being  used  by  the  MIT  led  consortium. 

The  model  construction  described  in  this  report  relies  on  hybrid  models  constructed  from  several  global  models 
of  the  crust  and  upper  mantle.  In  work  during  the  second  year,  University  of  Connecticut,  in  cooperation  with 
the  other  members  of  our  consortium,  starting  hybrid  models  will  be  refined  to  be  consistent  with  ground  truth 
travel  times  being  assembled. 
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Figure  1.  Three-dimensional  crust  and  upper  mantle  structure  beneath  Nilore,  Pakistan  (IMS  station  NIL) 

assembled  from  the  Cornell  Eurasian  Moho,  Crust  5. 1,  and  the  RUM  upper  mantle.  Note  the  tetrahedral 
parameterization  in  the  uppermost  4  km  needs  to  be  densified  to  better  describe  shallow  structure  and 
eliminate  artifacts. 
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Figure  2.  P  Ray  endpoints  in  the  vicinity  of  NIL  for  rays  shot  at  constant  increments  of  vertical  and  azimuthal 
take-off  angles.  The  Cornell  Eurasian  Moho  topography  is  superposed.  The  hybrid  3-D  model  in  which 
rays  were  shot  was  parameterized  by  tetrahedra  of  variable  size  and  its  size  is  25°  lat  x  25°  long  x  700 
km  deep.  Spacing  of  knots  at  tetrahedron  vertices  ranged  from  10  to  50  km. 
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